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Abstract 

Neurons are the central biological objects in understanding how the brain works. 
The famous Hodgkin-Huxley model, which describes how action potentials of a neuron 
are initiated and propagated, consists of four coupled nonlinear differential equations. 
Because these equations are difficult to deal with, there also exist several simplified 
models, of which many exhibit polynomial-like non-linearity. Examples of such models 
are the Fitzhugh-Nagumo (FHN) model, the Hindmarsh-Rose (HR) model, the Morris- 
Lecar (ML) model and the Izhikevich model. In this work, we first prescribe the 
biologically relevant parameter ranges for the FHN model and subsequently study the 
dynamical behaviour of coupled neurons on small networks of two or three nodes. To do 
this, we use a computational real algebraic geometry method called the Discriminant 
Variety (DV) method to perform the stability and bifurcation analysis of these small 
networks. A time series analysis of the FHN model can be found elsewhere in related 
work |15) . 

The Fitzhugh-Nagumo Model The Fitzhugh-Nagumo (FHN) system of equations is 
a prototypic model of excitability. Introduced by FitzHugh |13j with equivalent circuit by 
Nagumo et al.[25j the system is a generalization of the Van der Pol oscillator [TT] and a 
reduction of the neural electrophysiological model due to Hodgkin and Huxley (HH) [IB] . 
As a generic model of excitability and oscillatory dynamical behavior, the FHN system is of 
relevance for a range of physical and physiological research topics pT | [?fl [TOl [51 H^. the 
most extensive of which are those concerned with cardiac [20 , 6 , 28J and neuronal |H] 1301 133] 
cell dynamics. The corresponding equations for n electrically coupled FHN neurons are 

dx Xj ^ \ , . dy^ / , \ / 1 \ 

=2;;- -^-yi-|-g2^(a;, -Xj), — ^ e{x^ + a - by,) (1) 

for i = 1, . . . ,n and where we have taken a e [—2, 2], b £ (0, oo), ~l < g < I {g ^ 0) and < 
e < 0.1. It is well-known that in the FHN model, the variables have no direct physiological 
interpretation. However, for the parameter ranges quoted above, the qualitative behaviour 
of the x's and y's are similar to that of the voltage and gating variables in the Hodgkin- 
Huxley equations. 



t Correspondence to: IRCCyN, 1, rue de la Noe, BP 92101, 44321 Nantes Cedex 03 France. Tel : -|-33 2 
40 37 69 00. Fax : +33 2 40 37 69 30 



stationary points in the FHN system Neurons within the central nervous system 
can exist in a variety of dynamical states. Many, for example, are in a state of quiescence 
and elicit a relaxation oscillation, known as an action potential, when perturbed with a 
suprathreshold stimulus. This is oft termed as "spiking". With the appropriate choice of 
parameters, the FHN model can exhibit such behaviour. 

However, upon varying these parameters, one can also invoke a Hopf bifurcation resulting 
in relaxation oscillations at an intrinsic frequency without the need for any stimulation. 
Such self sustained neurons can be found within the central nervous system and can have 
complex interactions with the environment. A typical example is that of circadian cells 
which act like the organisms clock cells and have been shown to have the ability to entrain 
their oscillations with the environments light-dark cycle |31j . 

Neurons have also been shown to express bistability [25], which again is present in the 
FHN model due to the cubic nonlinearity present in the system. 

The dynamical states in the FHN model described above can be identified via a stability 
analysis of the stationary points in the system corresponding to ^ = iff" = The 
solutions to this system shall henceforth be referred to as simply the steady states. Such 
an analysis is thus essential if one is to have a basic understanding of the system. In fact, 
the steady state equations of this model have already been solved exactly for the n ~ 1 
and n — 2 cases Recently a conjecture relating the steady states of a system and 
its Kolmogorov-Sinai entropy has brought to greater hght the importance of steady state 
analysis for neural modeHng in general jH [TJ |51 1^. 

It should be noted that computational and numerical algebraic geometry methods have 
found many applications in many branches of theoretical physics in general (see, e.g., 
Refs. [231 1231 El])- Below we use the DV method to solve the corresponding equations 
for the n = 3 case and study the stability and bifurcation structure of these systems. The 
functions we use are in the packages Groebner and RootFinding [Parametric] of the com- 
puter algebra system Maple 13. Due to the polynomial-like nonlinearity also exhibited 
in the HR, ML and Izhikevich models, the same technique may also be applied to these 
systems. 

The problem addressed in this paper can also be related to the more general problem 
of the algebraic analysis of the solutions of a differential system, studied for example in |7], 

mm or mm- 

It is also worth mentioning that, though we have performed our computations in Maple, 
some of the computations (quantifier elimination, sarnple points extraction) can be run 
using Discoverei0 QEPCAeQ REDUCE^ Mathematicaly STRINGVACUjj^] or RA^ 

Algebraic tools For this study, we used the Discriminant Variety [TSJ and Cylindrical 
Algebraic Decomposition [TD]. The Discriminant Variety is an implicit representation of 
the desired partition, while Cylindrical Algebraic Decomposition describes explicitly each 
cell of the partition. 

Combined together, these two methods are well adapted to the analysis of the steady 
states (stable steady states). They provide a partition of the parameter space into con- 

*http://www. is. pku.edu.cn/ xbc/discoverer.html, developed by Wang, Xia et al. 
thttp://www.usna.edu/Users/cs/qepcad/B/QEPCAD.html, developed by Hong et al. 
thttp://reduce-algebra. sourceforge.net/, developed by Hearn, Weispfenning et al. 
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Figure 1: Parameter space 
decomposition adapted to 
the number of steady states 
for n ~ 3: in each 
connected component out- 
side the surface, there is a 
constant number of steady 
states. 
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Figure 2: Number of steady and stable states for a — 0, 
b — 2 according to the parameter g. 



nected cells, such that within each cell, the number of steady states (stable steady states) 
is constant. 

Finally, the main algebraic criterion to decide if a solution is stable is the Routh-Hurwitz 
criterion, or its Lienard-Chipart variant |17j . For our problem, we used a reduced criterion, 
more adapted to the algebraic DV and CAD methods. 

Results The model in Eq. [l] has already been studied for the case n = 2 in Ref. [9j . We 
therefore focussed our work on the case n > 3. 

First, we succeeded in describing the steady steady states for n = 3. The result of our 
computation is summarized in Figure [l] 

However, when n > 3, the computations of the stable states are much more difficult 
and we did not succeed in describing the full parameter space according to the number 
of stable steady states. Nevertheless, by fixing the parameters a and b to the numerical 
values (respectively and 2), we were able to describe the number of steady and stable 
steady states according to the free parameter g, for the n — 3 case. Since there is only one 
parameter, the description of the parameter space is a union of intervals in the parameter 
g for which the number of stable (and steady) solutions is constant. 

The results of the computations are summarized in Figure [2] 
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